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We consider self-dual transverse-field Ising spin chains with m-spin interaction, where the phase 

transition is of second and first order, for m < 3 and m > 3, respectively. We present a statistical 

analysis of the spectra of the Hamiltonians on relatively large L < 18 finite lattices. Outside the 

critical point we found level repulsion close to the Wigner distribution and the same rigidity as for 

the Gaussian Orthogonal Ensemble. At the transition point the level statistics in the self-dual sector 

is shown to be the superposition of two independent Wigner distributions. This is explained by the 

existence of an extra symmetry, which is connected to level crossing in the thermodynamic limit. 

Our study has given no evidence for the possible integrability of the models for m > 2, even at the 

transition point. 
PACS numbers: 05.50.+q, 05.30.-d, 05.20-y, 05.45+b 

I. INTRODUCTION 

There are many problems in physics in which multi-particle interactions play an important role. One may mention 
nuclear forces, solid 3 He Q, adsorbed systems [|| and plasmas ||. It is known from some exact results [Q that 
the critical properties of models with many body forces generally depend on the range of interaction. Recently 
considerable effort has been made to clarify the properties of a simple one-dimensional quantum model described by 
the Hamiltonian ||||: 

n = -jJ2 ofof +1 • . . af+m-x -hJlaf ^-JH X - hH z . (1) 



Here the erf, cr ; z are Pauli matrices at site I and J and h are the exchange coupling and the transverse field, respectively. 
The classical statistical mechanical equivalent of this model is a two-dimensional square lattice Ising model with mixed 
m-spin and two-spin interactions 0. 

The Hamiltonian Eq. ([!]) is self-dual @|| and the self-dual point is J = h independent of m. According to numerical 
studies [|u]-|l4|| there is one phase transition in the system, which takes place at the self-dual point, and the transition 
changes from second to first order, when m > m c = 3. In the borderline case m — 3 the transition is conjectured |7] 
to belong to the four state Potts universality class, a conjecture which is supported by an approximate mapping |1E | 



and by numerical studies M,[l^,|r4 16 



Concerning the simple structure of the model, its self-dual symmetry and the expected relation to Q-state Potts 
models, one can also pose the question, whether the model is integrable, at least in its self-dual point. To find an answer 
to this question in this paper we are going to study the statistical properties of the spectrum of the Hamiltonian. As 
$H ' it has been established in a series of papers fl7r]23]l the spectrum of a Hamiltonian (or the transfer matrix for classical 
statistical mechanical models) has different statistical properties for integrable and non-integrable models and one 
can make a close connection to the theory of the spectral properties of random matrices. In the actual calculation 
we first make use of all all those symmetries of the Hamiltonian which do not depend on the value of the couplings, 
and for large finite lattices we block-diagonalize the eigenvalue matrix of the problem. The statistical analysis of the 
energy levels is then performed for each block separately. In a non-integrable model, in which no further internal 
symmetry is present the matrix-elements of a block-matrix are expected to be loosely correlated, so that they can be 
approximately represented by random entries. Indeed the spectrum of non-integrable models is found to belong to 
the class of orthogonal random matrices, to the so called Gaussian Orthogonal Ensemble (GOE) and the level spacing 
distribution is described by the Wigner surmise |24[: 
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P( s ) = | S exp(-7T S 2 /4). (2) 

On the other hand if the Hamiltonian is integrable by the Bethe ansatz there is an infinite number of internal sym- 
metries and consequently the matrix-elements of a block-matrix are strongly correlated. Loosely speaking integrable 
Hamiltonians are so peculiar that they are not well described by an 'average Hamiltonian'. Then one expects that 
in this case the eigenvalues themselves behave like independent random numbers, so that the spectrum of integrable 
models belongs to the ensemble of diagonal random matrices and the level spacing distribution is described by the 



Poissonian (exponential) distribution: P(s) — exp(— s). Numerical studies of integrable models |17-23| are indeed in 
agreement with this assumption. 

In this paper we are going to perform the analysis of the level statistics of the multispin coupling Hamiltonian in 
Eq. ([!]). We are going to answer two questions. The first question is whether or not the Hamiltonian is integrable, at 
least at the transition (self-dual) point. Our second question concerns the characteristics of the level distribution at a 
first order transition point. The paper is organized as follows. The symmetries of the Hamiltonian in Eq. ([!]), which 
are essential to perform a block-diagonalization, are presented in Section 2. The statistical analysis of the spectrum 
of the block-diagonalized Hamiltonian is given in Section 3, while the results are discussed in the final Section. 



II. SYMMETRIES OF THE HAMILTONIAN 



As described in the Introduction the first step in a statistical analysis of the energy eigenvalues is to block-diagonalize 
the Hamiltonian in Eq. (|l|) using all those symmetries of the model which do not depend on the actual values of the 
couplings. Before analyzing these symmetries, let us first notice that if E{J, h) = {Eo(J, h), E\ (J, h), ■ ■ ■ , E 2 l_i(J, h)} 
denotes the set of energies of the Hamiltonian Eq. (|l]), one has E{J, h) = E{±J, ±h). This can be seen introducing the 

operators O x — Y\^I^ of and O z = [j^™ -1 af m and noting that H a O^ — eapO^Ha, where a,/3 = x,z and e Qj/ g = 1 
if a = /3 and -1 otherwise. On a finite lattice this symmetry holds for periodic boundary conditions and if the length 
of the chain is a multiple of to. In what follows we consider this type of lattices and restrict ourself to the case h > 
and J > 0. 

The symmetries of the model are of three types: i) space-like symmetries, which describe invariance of the system 
under geometrical transformations (translation, inversion, etc) ; ii) gauge symmetries, which are connected to invariance 
of the Hamiltonian under internal transformations; and finally iii) duality symmetries, which make a connection 
between the strong- and weak-coupling regimes of the Hamiltonian. 

i) The space symmetry of the model on a finite lattice depends on the boundary condition. In a statistical analysis 
of the spectrum of finite systems it is desired to use the most symmetric boundary condition to get a block structure, 
which well represents the statistical behavior of the spectrum in the thermodynamic limit. Therefore, as already 
mentioned, we apply periodic boundary conditions, which can be formally expressed as er£ +i = erf- The space 
symmetry group is then the automorphy group of a ring, irrespective of the range of the interaction to. This is the 
dihedral group generated by the translation T and the reflection R, (T L — R 2 = Identity and TR — i?T L_1 ) which 
both commute with TL. 

ii) The gauge symmetries are generalizations of the spin-reversal symmetry for the well known to = 2 case. Recalling 
that we take L to be a multiple of to, let us introduce a set of n = 2 m ~ 1 operators Ok for k = 0, 1, . . . , n — 1: 

L/m—l am+m — 1 

o k = n n ( 3 ) 

a— i—am 

where ho, ki, ■ • • , &m-2 are the bits of the binary representation of k, and fc m _i is such that 

rn 1 

ki even . (4) 

i=0 

The spin reversal symmetry, which holds when m is even, is O n _i corresponding to ki — 1 for all i. It is straightforward 
to check that the condition Eq. ^ ensures that all the operators Oi commute with TL{J,h). It is also clear that these 
operators are diagonal, involutive and form an Abelian group (Oq is the identity). This implies that all the 2 m ~ 1 = n 
representations are one-dimensional and the corresponding projectors are of the form Pr = 1/n y^Lq efO^ where 
ef = ±1 for all i and R. All these projectors split the Hilbert space in 2 m_1 invariant subspaces of size 2 L ~ m+1 each. 
For example, for to = 2 and L even, Pq — (l/2)(Oo + 0i) projects onto the subspace with an even number of up spins, 
while Pi = (1/2)(Oq — Ox) projects onto the subspace with an odd number of up spins. The projector Pq = l/n^Oi 
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projects onto the most symmetric subspace to which the ground state belongs (we refer now to this subspace as the 
ground-state sector), whereas the other 2 m ~ 1 — 1 sectors become degenerate in the thermodynamic limit. Thus, in 
this limit, the degeneracy of the ground state in the strong coupling phase J > h is given by 2 m ~~ 1 . This degeneracy 
for the m = 3 model is just four, which led Debierre and Turban M to conjecture the same universality class for the 
transition as that for the Q — 4 state Potts model. 

The combination of the space symmetry and the gauge symmetry is not obvious, since the operators of these two 
groups do not commute in general. The product of these two groups is a semi-direct product, (not a direct product) 
since the gauge group is a normal subgroup. As usual the states are labelled by the number of the representation, R, 
to which they belong. We have computed the character table of the symmetry group from which the dimensions of 
the invariant subspaces are deduced and then the block-diagonal Hamiltonian is constructed. The dimensions of the 
irreducible representations and the size of the corresponding blocks are given in Tables I-III for different values of m, 
in the range of < R < L/2 + 3 and < R < (L — l)/2 + 2, for L even and odd, respectively. We note that in the 
ground state sector, which is labelled by R — and corresponds to Po = (1/n) J^Oj all the operators of the space 
and gauge symmetry group commute, so that in this sector we have a representation of the dihedral group. In what 
follows we use the same labelling convention as in J|^] . 

iii) As mentioned in the Introduction the Hamiltonian in Eq. ([lj) has the property of duality symmetry. To show 
this and its consequences in finite lattices, first we define, for an infinite lattice, dual Pauli operators rf , rf as 

_ Z ^X X X / r\ 

Z XX X (r>\ 

a i — ~i ~i+l ■ ■ ■ T i+m-l \°) 

in terms of which the Hamiltonian in Eq. (Q) is expressed as 

n = J E T i h Y. T i T i+i ■ ■ ■ T i +m -i • (7) 



Consequently the two sets of energies E{ J, h) and E(h, J) are equal: 

E(J,h) = E{h,J) (8) 

and the self-dual point h = J corresponds to the transition point of the system, provided there is one single phase 
transition in the system. The duality symmetry, as described above holds in the thermodynamic limit, i.e. when the 
length of the system L — > oo. In a finite system duality generally relates sectors of the Hamiltonian with different 
boundary conditions. With periodic boundary conditions one has the symmetries er£ +i = erf and erf, 1 = af, which 
in terms of the dual operators in Eq. (||) and Eq. are only satisfied in the ground state sector of the Hamiltonian. 
As a result self-duality holds only in the ground state sector, which is indeed verified numerically. Based on this 
observation we expect somewhat different statistical properties of the energy levels in the self-dual and non-selfdual 
sectors. 



III. RESULTS OF THE RANDOM MATRIX THEORY 



Using the symmetries as described in the previous Section we have performed the block-diagonalization of the 
eigenvalue matrices for large finite lattices, the size of which was a multiple of the length of the interaction m. We 
went up to L = 18, 16 and 15 for m = 3, m — 4 and m — 5, respectively. The size of the blocks, as seen in Tables 
I-III, is relatively small, especially for larger values of m their size is reduced by gauge symmetry. 

Having the block-diagonalized Hamiltonian we solved their spectrum by standard numerical methods, which are 
contained in the LAPACK library. The next step, before performing the analysis, is to unfold the spectrum, i.e. to 
subtract the average tendency and to keep only the fluctuations, which are normalized in the same manner at each 
part of the spectrum. Technical details relating to unfolding the spectrum are given in Ref ^||,^3| . 

The unfolded spectrum is then analyzed and several spectral quantities are determined and compared with the 
predictions of random matrix theory. First, we consider the level spacing distribution, P(s), which is expected of the 
Wigner form in Eq. (|^) for non-integrable models, whereas it is generally of the Poissonian form for integrable models. 
To analyze realistic spectra it is often useful to consider Brody's interpolation formula: 



Pp(s) = c(l + /3)s^exp (-cs^ 1 ) 



(9) 



with c 



r 



(§+2\ 



1+/3 



, which corresponds to the Wigner and the Poisson form for f3 = 1 and (3 = 0, respectively. 

The interpolation parameter (3, which is determined by an optimalization fit, proved itself to be a useful indicator for 
the localization of integrable varieties 22 23| ]. 
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Another quantity characterizing the independence of the eigenvalues is the spectral rigidity in an interval of length 



where N u (e) = J^i ®( e — e i) is the integrated density of unfolded eigenvalues and (. . .) a denotes an average over a. 
Finally, we shall also consider the number variance £ 2 (Z) defined as the variance of the number of unfolded eigenvalues 
in an interval of length I: 



where the brackets denote an averaging over e. 

FIG. 1. Level spacing distribution for L=15 M=3 h/J=1.36 and ALL the representations. The exponential ("E") and the 
Wigner ("W") distributions are shown (full line), together with the Brody distribution (broken line) for the fitted best value 
of the parameter /3 = 1.01 (see text). 

First we presents the results of the statistical analysis of the spectra outside the transition point. As seen in Figs, [j] 
and U on the example of the m — 3 model at a coupling J/h — 1.36 for a 15-site chain, all the three characteristic 
quantities of the spectrum are very well described by the Wigner distribution p4}] . Fig. || presents the rigidity and 
the number variance for the same parameters. The expected behavior for independent random energies and for the 
eigenvalues of the GOE are also shown. The GOE behavior is observed up to quite large values of L, indicating that 
GOE matrices provide a good description of the Hamiltonian. The data shown are obtained averaging over all the 
representations. However very similar results are obtained averaging only over the self-dual sector. The characteristic 
parameters of the spectrum do not depend on whether the sector under consideration is selfdual or non-selfdual We 
also note that very similar behavior is found for other ranges of the interaction m > 3 or for other values of J/h > 1. 
Consequently there is no evidence for the integrability of the model with m > 2 away from the critical point. 

FIG. 2. Rigidity (a) and variance (b) for L=15 M=3 h/J=1.36 and ALL the representations. The corresponding quantities 
for the GOE matrices and for random diagonal matrices are also shown in full lines. 

In the following we investigate the level statistics of the model as a function of the ratio h/J and calculate the 
interpolation parameter (3 in Eq. (||) as a best fit over the self-dual and non-selfdual sectors. The results are shown 
in Fig. H for L = 15 and m = 3, whereas data for the largest system size are only included in the self-dual sector. 
We note that the corresponding data for m = 4, m = 5, and a less extensive calculation for m = 6, lead us to very 
similar conclusions. One can see in Fig. || that Brody's parameter (3 has different behavior in the self-dual and in the 
non-selfdual sectors. While in the non-selfdual sectors f3 is approximately constant and its value (3 w 1 corresponds 
to the GOE result, in the self-dual sector there is a change in the value of (3 around the self-dual point. Actually its 
value drops from (3 — 1 to about (3 « 0.45 at the transition point. The region, where the change in (3 takes place 
seems to shrink only to the self-dual point in the thermodynamic limit, as can be seen in Fig. by comparing the 
results with L = 15 and L = 18. 

FIG. 3. Parameter /3 as a function of h/J for m = 3. For L = 15 the data are averaged separately over all self-dual sectors or 
over all non self-dual sectors. For L — 18 only the representation Rq to which the ground state belongs is taken into account. 

This observation leads us to study carefully the spectral properties of the models at the transition point, the results 
of which are shown in Fig. |](a). As seen in the figure the level spacing distribution could not be well fitted by 
the interpolation formula in Eq. (j^), at least with the symmetries we have taken into account. However, it is given 
approximately by the arithmetic average of the Wigner and Poisson distributions, which is also shown in the figure. We 
argue that the measured spectral quantities in the self-dual sector can be interpreted as if the spectrum is composed 
of two independent Wigner distributions. To check our assumption we have taken two non-selfdual blocks of roughly 
the same size each of which has Wigner characteristics and merged the levels of the two blocks. Then we analyzed 
the level statistics of the combined blocks and the obtained results in Fig. [|(b) looks very similar to those we found 
for the self-dual sector and presented in Fig. £|(a). 

FIG. 4. a) Level spacing distribution for L=18 M=3 and for the representation to which the ground-sate belongs, 
b) Combination of the spectrum of two different representation Ri$ and Ris of the m = 5 model for L — 15. 



/: 




(10) 




(11) 



4 



Thus at this point we conclude that the spectrum of the self-dual sector at the self-dual point is seemingly composed 
of two independent parts, each having Wigner-type characteristics. This type of behavior is the result of an extra 
symmetry, the self-duality, which is just seen in the selfdual sector. Furthermore, we argue that this extra symmetry 
is manifested by the crossing of energy levels at the self-dual point in the thermodynamic limit L — > oo. To see this 
we have calculated the quantity: 

\ nun(sj, Si-i) I i 

which measures the asymmetry in the the level spacing distribution. For independent random variables chosen 
according to a Wigner distribution one has: P{5) — rrjfzp yielding a mean value < 5 >= 1 + n/2 w 2.5708. For 

matrices from the GOE, the correlation between spacings Si modifies this value. We have numerically found [] that 
this value does not vary considerably with the size of the matrix, and is very close to 

< 5 >goe~ 3 (13) 

As seen in Fig. | for non-selfdual sectors 5 is indeed close to the GOE result in (|l3|). Similarly, S w 3 is found in the 
self-dual sector far from the self-dual point, however there is a sharp increase in S in the neighborhood of h/J — 1. 



[ |25[ Since the value of S at the self-dual point is monotonically increasing with the size of the system (see Tab. IV), 
one expects that in the thermodynamic limit 5 — ► oo. Thus there is an exact degeneracy in the selfdual sector at the 
transition point, which should be connected with the presence of an extra symmetry. The possible origin of this extra 
symmetry is discussed in the final Section. 

FIG. 5. The asymmetry parameter 5 in Eq. ( p^ ) as a function of h/J for L = 16 and m = 4. The upper curve corresponds 
to the self-dual sector and the lower curve to non self-dual sectors. The error bars are calculated as the variance divided by 
the square root of the number of spacing ratios. 



IV. LEVEL STATISTICS AT A FIRST-ORDER TRANSITION POINT 



In this paper we have studied the statistical properties of the spectrum of a transverse-field Ising spin chain with 
771-spin interactions and compared to the predictions of random matrix theory. Away from the transition point, 
which is known exactly from duality symmetry, the spectrum is shown to be a Gaussian Orthogonal Ensemble and 
its properties are well described by the Wigner distribution. On the other hand at the transition point the spectral 
properties of the selfdual and non-selfdual sectors are different. While the spectra of non-selfdual sectors are close 
to the Wigner distribution the same for the selfdual sector can be described as the composition of two independent 
Wigner distributions. Furthermore, we have shown that this special behavior at the transition point is the result of 
level crossing in the thermodynamic limit. 

This observation can be compared with the known exact |26| and numerical |2?J results on the two-dimensional 
Q > 4 Potts model. As known exactly ps[| , this model is also self-dual and there is a first order transition in the 
system. As an analogous quantity to the Hamiltonian in (Eq. [I]) we consider the T transfer matrix of the Potts model, 
which in the Hamiltonian limit ]29| , |3"c| is given by T — exp(— rHp), where r denotes the lattice spacing and Hp is 
the Hamiltonian of the one-dimensional quantum Potts model. According to exact results p6| in the thermodynamic 
limit the ground state of Hp at the transition point is (Q + l)-fold degenerate. At this point the first two levels of 
the self-dual sector, as well as the first levels of the Q — 1 other, non- selfdual sectors are degenerate. Thus the first 
order nature of the transition is manifested by a level crossing in the self-dual sector. (For finite rings one observes 
a hybridization gap in the selfdual sector, which vanishes exponentially with the size of the system |27| .) As shown 
by numerical calculations the same type of level crossing phenomena takes place for the higher lying levels, too. 
Thus, for finite systems, the spectrum at the transition point is expected to decompose into two parts, which are 
going to be degenerate in the thermodynamic limit. 

Our numerical results on the multispin coupling model are in agreement with the above picture, thus we expect a 
similar scenario. The selfdual symmetry at the transition point, which is connected to a level crossing in the selfdual 
sector in the thermodynamic limit is responsible for the unusual spectral properties of the multispin coupling models 
for to > 3. The m = 3 model, in which the transition is expected to be second order, is assumed to represent the 



Miagonalizing 1000 GOE random matrices of size ranging from 3 up to 2000. 
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border limit of continuous models. Thus one expects that the above scenario, which stays valid as m — > 3 + , could 
hold also for m = 3, perhaps with another type of functional form for the size dependence of the hybridization gap. 
Indeed, this type of behavior has been found by our numerical studies. 

Finally, we turn to discuss possible integrability of the multispin coupling model in Eq. (Q). As a result of our 
numerical studies of spectral properties of the model we conclude that there is no evidence in favor of integrability of 
the Hamiltonian in (Eq. |l|) for m > 2, even at the transition point. 
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TABLE I. Dimensions of the irreducible representations (i.e. degeneracy) and size of the corresponding block. 



M=4 L=16 


label R 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


16 


dimension 


1 


1 


1 


1 


1 


1 


1 


1 


2 


2 


2 


2 


2 


2 


2 


2 


2 


size 


330 


265 


202 


265 


288 


288 


224 


224 


529 


512 


544 


480 


544 


480 


496 


512 


526 


label R 


17 


18 


19 


20 


21 


22 


23 


24 


25 


26 


27 


28 


29 


30 


31 


32 


33 


dimension 


2 


2 


2 


2 


2 


2 


2 


2 


2 


4 


4 


4 


4 


4 


4 


4 


8 


size 


512 


496 


512 


496 


512 


526 


512 


496 


512 


1088 


960 


1024 


1024 


1024 


1024 


1024 


2048 



TABLE II. Dimensions of the irreducible representations (i.e. degeneracy) and size of the corresponding block. 
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TABLE III. Dimensions of the irreducible representations (i.e. degeneracy) and size of the corresponding block. 





L=10 


L=12 


L=15 


L=16 


L=18 


m=3 




4.79 


6.54 (9.63) 




11.31 


m=4 




4.58 




5.09 (8.67) 




nr=5 


1.76 




4.33 







TABLE IV. The asymmetry parameter S in Eq.(|l2j) as a function of L for m = 3, 4 and 5 in the ground-state sector. When 
present the number in parenthesis refers to the asymmetry calculated for the entire self-dual sector. 
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